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Abstract. A reduced Keller-Segel equation (RKSE) is a parabolic-elliptic system of 
partial differential equations which describes bacterial aggregation and the collapse of 
a self-gravitating gas of brownian particles. We consider RKSE in two dimensions, 
where solution has a critical collapse (blow-up) if the total number of bacteria exceeds 
a critical value. We study the self-similar solutions of RKSE near the blow-up point. 
Near the collapse time, t = t c , the critical collapse is characterized by the L oc (tc — t) 1 / 2 
scaling law with logarithmic modification, where L is the spatial width of collapsing 
solution. We develop an asymptotic perturbation theory for these modifications 
and show that the resulting scaling agrees well with numerical simulations. The 
quantitative comparison of the theory and simulations requires to take into account 
several terms of the perturbation series. 
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1. Introduction 

In this paper we consider a reduced Keller-Segel equation (RKSE) 

dt P = Ap-V-(pVc), 
Ac = -p, 

which is the parabolic-elliptic system of partial differential equations for two scalar 
functions, p = p(r,t) and c = c(r, t). Here r £ f2 C IR- is the spatial coordinate in 
dimension Z) and £ is the time. We assume that either Q = IR D or Q is a bounded 
domain. For Q = IR D , we also assume that both p and c decay to zero as |r| — > oo. In 
the bounded domain case, we assume the zero flux condition for both p and c through 
the boundary dfl. 

RKSE is the reduction of the well-known Keller-Segel model (also sometimes called 
Patlak-Keller-Segel model). See e.g. [HEIOIIIIEIIHIITIIHIIHI HOI HH H21 HH H5l US] 

and references therein. The Keller-Segel model was derived for the macroscopically 
averaged dynamics of bacteria and biological cells. Below we refer to bacteria and cell as 
synonyms. Bacteria often communicate through chemotaxis, when bacteria both secrete 
a substance called chemoattractant and move along the gradient of chemoattractant. 
The macroscopically averaged dynamics of bacteria is described by the bacterial density 
p(r,t) and the chemoattractant concentration c(r, t). Bacteria are self-propelled and, 
without the chemotactic clue, the center of mass of each bacteria typically experiences 
a random walk. The random walk is described by the first term (diffusion) on the 
right-hand side (rhs) of the first equation in (CQ). The diffusion of chemoattractant 
is described by the Laplacian term in the second equation. The term on the rhs of 
the second equation corresponds to the production rate of chemoattractant by the 
bacteria, which is proportional to the bacterial density. The second term on the rhs 
of the first equation characterizes the motion of bacteria towards large values of c. The 
motion of bacterial colonies is thus determined by competition between random-walk- 
based diffusion and chemotaxis-based attraction. For the convenience of the readers, 
we provide a more extensive description of the Keller-Segel model and the derivation of 



RKSE in Appendix A 



Equation (pD) also describes the dynamics of a gas of self-gravitating Brownian 
particles, which has applications in astrophysics including the problem of stellar 
collapse [T71 HH HH [19]. In this case, the second equation in (JTJ is the Poisson 
equation for the gravitational potential, — c, while p is the gas density. (All units are 
dimensionless) . The first equation in ([1]) is a Smoluchowski equation for p. Below we 
refer to p and c as the density of bacteria and the concentration of chemoattractant, 
respectively, but all results below are equally true for the gravitational collapse of a gas 
of self-gravitating Brownian particles. 

A solution of RKSE in dimension one (D=l) is global (in time). For D > 2 (e.g. 
in dimensions two and three ) a finite time singularity occurs [6] provided the initial 
condition is large enough. E.g. for dimension two (D=2), a finite time singularity occurs 
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for N > 8tt, where N = J p dr is the total number of bacteria (in rescaled units) 
Below we focus on unbounded domains Q = M. D . 

The formation of singularity in a finite time (blow up) is a quite general phenomenon 
observed in many nonlinear systems including self-focusing in nonlinear optics, plasmas, 
hydrodynamics, and collapse of Bose-Einstein condensate [21]. Blow up is often 
accompanied by a dramatic contraction of the spatial extent of solution, which is then 
called by collapse [21]. Collapse typically occurs when there are (i) self-attraction in 
nonlinear systems and, (ii) a conserved quantity, such as the spatial norm (e.g., L2 or L\ 
norm) of the solution. Such systems are often described by the nonlinear Schrodinger 
equation (NLSE) [22]: 

id t ij + V 2 ^+ |^| 2 ^ = 0. (2) 

for complex variable ijj(r,t). NLSE conserves the integral P = J \tp\ 2 dr and supports 
collapse for D > 2. Similarly, RKSE describes the attraction between brownian particles 
and conserves L\ norm of p as well as RKSE admits collapse for D > 2. 

A collapse in RKSE corresponds to the aggregation of bacterial colonies in 
biological applications and gravitational collapse for self-gravitating Brownian particles. 
Aggregation is a first step to a formation of multicellular organisms and quite important 
in biological applications [TJ. E.g., the evolution of a low-density Escherichia coli 
bacteria colony in a petri dish is about one day [5]. However, if the bacterial density 
is locally high then bacteria aggregate on a timescale of several minutes [5]. Thus 
the aggregation has an explosive character (see more details on that in Appendix 



Appendix A). Near singularity the Keller-Segel model is not applicable when typical 



distance between bacteria is about or below the size of bacteria. In that regime a 
modification of the Keller-Segel model was derived from microscopic stochastic dynamics 
of bacteria [231 EH EH]. That modified model prevents collapse due to excluded volume 
constraint (different bacteria cannot occupy the same volume). Here however the 
original RKSE without regularization is considered. 

Collapses in NLSE and RKSE have much common, as detailed in Ref. [H]. E.g., 
the number of particles P in NLSE has a similar meaning to the number of bacteria N in 
RKSE. One can also recall that \ip\ 2 is the probability density in quantum mechanics. In 
two dimensions (D = 2), the critical number of particles, P c = 11.70 .. . (for NLSE), or 
the critical number of bacteria, N c = 87r (for RKSE), determine the boundary between 
collapsing and noncollapsing regimes in both systems [61 [261 EH EH E3 [301 EI] • Collapse 
in the critical dimension D = 2 is strong for both RKSE and NLSE, which means that 
a finite number of bacteria (particles) is trapped within the collapsing spatial region. 
For the supercritical case (D > 2), collapse in both RKSE and NLSE is weak. Weak 
collapse implies that the collapse is so fast that particles (bacteria) cannot keep up with 
the collapse rate. Then a vanishing number of bacteria (particles) are trapped inside 
the collapsing region in the limit t — > t c . 
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1.1. Summary of results 

In this paper, we focus on the 2D self-similar solution of RKSE, Eq. (JI)). We assume 
that the spatial location of the collapse is r = 0. Near t c , in the neighborhood of the 
collapse, the solution has the following radially symmetric form: 

1 8 

P 



L(t)2(l + y 2 ) 2 
c=-21n(l + 2 / 2 ), 



(3) 

r := r, 



* L(t) 
L(t) for t -> t c . 

Here, L(t) is the time- dependent spatial width of solution. We also refer to L(t) as 
the collapse width. (We sometimes omit the argument of L for brevity.) The self- 
similar form ((3]) is valid in the limit t — > t c in the small spatial neighborhood of the 
collapse point. This local applicability of the self-similar solution is typical for collapses 
in numerous nonlinear systems [21] . 

A number of different scalings for L(t) have been proposed. First is the scaling, 



L {t) = Cv ^-rt e -y Z ^P r [_l n (t c _t)](l/4)(-ln(* c -*))-- ; (4) 

where c is an unknown constant. This scaling was derived in Ref. [I] using formal 
matched asymptotic expansion of RKSE near fl3]). 
The second scaling, 



1/2 



2+-y , / \n(tc— t) 

Lit) = 2e~^ 2 Vt7 zr te"y ~, (5) 

was derived in Refs. [§] and [TJ]. Here 7 = 0.577216 ... is the Euler constant. In Ref. [9], 
the formal matched asymptotic expansion of RKSE was used. The approach in Ref. [T4~| 
was based on the expansion of the perturbation around the collapsing solution in 
terms of the eigenf unctions of the linearization operator. Refs. [9J and [H] give different 
estimates of errors. 
The third scaling, 

1 / ln(t c -t) In [- In (t c ^tj] 

L(t) = cy/t~^te~^y 3 , (6) 



where c is an unknown constant, was obtained in Ref. [TTJ by somewhat heuristic 
arguments. (Also see Ref. [32] for more discussion.) 

The scaling laws dlJ-Q share two main features: (i) the leading order square- 
root dependence, L[t) oc y/t c — t, and (ii) the logarithmic-type modifications of L(t). 
These modifications are necessary for the building the theory of the collapse in the 
critical dimension (2D). Both features are strikingly similar to the critical collapse in 
the 2D NLSE [261 12TI [28l EH [331 EH EH [361 EZl [38]. Earlier simulations [SJ E] show 
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that corrections to the leading order scaling L(t) oc y/t c — t are necessary, but fail to 
determine the form of the corrections. 

In this paper we go much beyond the accuracy of scaling laws (J4j) — (JSJ) - We derive a 
new scaling law that agrees with the direct simulations of RKSE. There are three main 
results in this paper. 

Our first main result is that L(t) is determined by the solution of the following 
ordinary differential equation (ODE): 

d T a 2 M b 



(I 



2 



ln± (ln±) 2 (ln±) 3 V(ln±) 4 



M = —2 — 27 + 2 In 2, (7) 

71-2 

6 = 21n 2 2 + 41n2 + 7(-4-27 + 41n2). 

3 

The adiabatically slow quantity 

a = -L(t)d t L(t), (8) 
evolves over a new time scale described by a new variable, r, defined as 

dt' 



T 



L(t>) 



(9) 



Here and below, the notation f(x) = 0(x) means that there exists a positive constant 
c such that |/| < c\x\ as x — > 0. It follows from OH]) that r — > 00 as t — > t c , so that 
r(t) maps the collapse time t = t c into r = 00 in full analogy with the "lens transform" 
of NLSE [321 HDJ EZ]- The decrease of L — > as t — > t c implies that a > 0, and the 
logarithmic modification of L(t) oc yjt c — t scaling results in a — > as t — > t c . The 
logarithmic modification also makes a a slow function of (t c — t) 1 ^ 2 , compared with L. 
These scalings, as well as the definition of a, are in qualitative analogy with the scaling 
for NLSE collapse. 

Our second main result is that the asymptotic solution of (JTj) in the limit t —> t c , 
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together with (jSJ) and OH]), is given by 

ln/3(t c -t) -l + &Inx 



2+7 f 

L(t) = 2e ~A/t c - texp i 



2 2x 



-1 + 2b + 2M(1 - ftlnx) , ( 1 \ + f( lnx 



,2 



x = y/-2ki P(t c -t) - M, 
M = -2 -7 + In 2, 

6=1 + ^, (10) 
o 



/3 = 2 exp < 2/* 




l l2 M+l 
/ = — In L — - In a H — In a - 

L := L(t ), a := a(t ) = -L(t )d t L(t ) 

This scaling was presented without derivation in Ref. jH]. The time t = t < t c 
is chosen arbitrarily, provided that at t — to, the solution is close to the self-similar 
form 0. (More details about choice of t are given at the end of Section \5\ and in 
Figured]) It is seen from (fTUj) that L(t) depends on the initial values L(t Q ) and d t L(to). 
The order of error terms in (TTOT) are discussed below, after Eq. ( 19T|) . 

Our third main result is the comparison of ( TTOT) with direct numerical simulations 
of RKSE. Figure [U shows excellent agreement between the theory and simulations. In 
the limit t — > t c , the new scaling f TTOT) reduces to (jSJ). We demonstrate, however, that 
while (jSJ) is asymptotically correct, it is in quantitative agreement with both (TIT)]) and 
simulations only for unrealistically small values 

L<1(T 10000 . (11) 

In contrast, the scaling f fTOj) is accurate starting from a moderate decrease of L(t) from 
the initial value L(0). Figure [2] shows the simulation with N = 1.0250A^ C , where (fit)]) is 
accurate (with the relative error < 7%) for L(t)/L(0) < 0.15. 



1.^. Outline of the paper 

The paper is organized as follows. In Section [2] we consider general properties of 
collapses in RKSE and their analogies with the collapses in NLSE. In Section [3] we 
study a collapsing self-similar solution of RKSE. We write the self-similar solution as 
a rescaled steady state solution in new "blow up" variables. In these variables, the 
self-similar solution transforms into the approximate steady-state solution. The full 
collapsing solution evolves slowly about the steady-state solution, and depends on the 
small adiabatically slow parameter a defined in (jSJ). We use a gauge transformation to 
a new dependent variable for perturbations about the self-similar solution. The gauge 
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Figure 1. Dependence L(t) obtained from the numerical simulations of RKSE (solid 
lines) is compared to the scaling (JSJ (dotted line) and to the scaling (JTUJ) (dashed-dotted 
lines). The lines of different colors correspond to different initial conditions (different 
values of N) . Different panels show the different orders of the scaling in the exponent 
of the first equation of (|10|) : (a) the terms up to O(x ) are taken into account; (b) the 
terms up to 0{x~ l ) are taken into account; (c) the terms up to 0(x~ 2 ), i.e., all terms 
except the error term 0(. . .), are taken into account. Convergence of the analytical 
results to the numerical results with increase of the order in inverse power of x is clearly 
seen in (a)-(c). The relative difference between numerical and analytical results in (c) 
is < 5% and decreases with the decrease of (N — N c ) /N c > 0. In simulations the initial 
conditions in the spatial Gaussian form as described in Section [7] 
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Figure 2. Dependence L(t) during the time interval significantly exceeding the time 
interval of self-similar regime. The solid line shows the results of the numerical 
simulations for N = 1.0250^, the dotted line shows the scaling ([5]), while the dashed 
line shows the scaling (fTU|) with all terms up to 0{x~ 2 ). The six-fold decrease of L 
from the initial value L(0) = 0.98773 already gives a good agreement between numerical 
simulation and (fT0|). with relative difference between them < 7% for L < 0.15. The 
scaling ([5]) agrees with simulation only in order of magnitude for L ~ 0.15. Figure [Th 
shows the same curves for N = 1.0250iV c zoomed-in to the origin. 



transformation brings the linearization operator about the self-similar solution to a self- 
adjoint form. In SectionH]we discuss the spectrum and eigenf unctions of the linerization 
operator. In Section [5] we expand the perturbations about the self-similar solution 
into eigenf unctions of the self-adjoint linearization operator, in order to derive a set 
of amplitude equations for the coefficients of the expansion. Compatibility conditions 
to ensure an adiabatic form of the expansion result in the ODE (J7|). In Section El 
we solve Eq. ((7j) to derive the scaling (flOl) . In Section [TJ we describe the simulation 
algorithm and the procedure for the extraction of the parameters of collapsing solutions 
from simulations. In Section IHJ the main results of the paper and future directions are 



discussed. In Appendix A we provide an extensive description of the Keller-Segel model 



and derive the reduced Keller-Segel equation. In Appendix B we provide the explicit 
expressions for the calculation of the scalar products from Section |5j these expressions 
are obtained by using the asymptotic expansions of the Meijer G-function and the F- 
function. 



2. Collapse of RKSE and NLSE 

Equation (CQ) has a form of a conservation law 

= - v ■ r, 



(12) 
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where V is the flux of the bacterial density given by 
T = — pV In p — 



(13) 



and c(r) is determined by the fundamental solution E(r, r') of the Poisson equation. 
The 2D case considered here implies that 



c r 



/ E(r, r'Wr'W, E(r, r) = — In |r - r'l 

J 27T 



(14) 



Eq. (fl4j) allows to rewrite Eq. ([I]) as a closed integro-differential equation for p. The 
integral term in that equation originates from ( fT4l) and represents the nonlocality of 
interaction due to diffusion of chemoattractant. 

Assuming decaying boundary conditions at infinity, we obtain the conservation of 
the total number N of bacteria: 

(15) 



N = J p(r)dr = const. 
One can also define a Lyapunov functional 



£ 



p(r) lnp(r) - p(r) 



p(r)c(r) 



dr. 



and represent Eq. in a gradient form 



M 



5£ 



dtp = V • ( pV j^J , — = \np-c, 

where the Lyapunov functional £ is a non-increasing function of time 
d£ 



(16) 



(17) 



dt 



— dr. 

P 



Functional £ is conserved only for a steady state solutions with zero flux T = 0. 

Although Eq. ([1]) is a gradient non-Hamiltonian system (as follows from ( fTTl) ). it 
has many striking similarities with NLSE ([2]) which can be written in a Hamiltonian 
form idtip 



4ir with the Hamiltonian 

dtp* 



H 



dr. 



(19) 



To prove existence of collapse in RKSE one can use a positive-definite quantity A 
J r 2 pdr, which determines the mean square width of bacterial density distribution 
4"2] . Vanishing of A guarantees the existence of collapse because of conservation of N. 
The proof of collapse existence for NLSE in D = 2 is based on a virial identity [27, 28]: 



8H. 



(20) 



Here B := J r 2 \ip\ 2 dr and H is defined in ( [19"]) . If H < then the positive-definite B 
turns negative in a finite time as follows from (|20l . It means that the negative value 
of the Hamiltonian is the sufficient condition for the collapse in NLSE. We also recall 
that in the quantum mechanical interpretation of NLSE, |-?/>| 2 is the probability density 
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of number of particles, i.e., the analog of p in RKSE. Thus A from RKSE is the analog 
of B in NLSE. However, RKSE is the non-Hamiltonian system and the direct analogy 
with a virial theorem for B does not work. Instead one can calculate a time derivative 
of A using Eqs. ([T|) and (fT4|) . integration by parts, and vanishing boundary conditions 
at infinity. For D = 2, this procedure gives: 



Here we also used symmetrization over r and r'. One concludes from (12 ip that A t < 
for N > 8tt and A turns negative in a finite time. That condition defines the critical 
number of bacteria 



because A < proves the existence of collapse by contradiction (A is the positive- 
definite). 

The existence of the critical number of bacteria (1221 is another similarity with 
NLSE, where the critical number of particles P c = J \tp\ 2 dr ~ 11.70 .. . The difference 
between collapses in NLSE and RKSE is that according to ( J2TT) for RKSE, any initial 
condition with N > N c develops into the collapsing solution in a finite time, while for 
NLSE, P > P c is the necessary condition for collapse but not the sufficient condition. 
Another qualitative difference between RKSE and NLSE is that RKSE is the integro- 
differential equation while NLSE is a partial differential equation (PDE). However, it 
was shown in Refs. [13|H3] that the generalized virial identity allows to prove the collapse 
existence in an integro-differential equation of NLSE-type with nonlocal nonlinearity. 
This type of nonlinearity describes, e.g., Bose- Einstein condensate with nonlocal dipole- 
dipole interaction. Collapse of such condensate was recently achieved in experiment jl5] . 

Qualitative similarities between collapses in RKSE and in NLSE can be also 
understood if we recall that RKSE is the mean-field approximation for the dynamics of 
self-gravitating gas of brownian particles, while NLSE is the mean-field approximation 
for the quantum dynamics of atoms with Bose statistics and attraction at ultra-cold 
temperatures. Thus both RKSE and NLSE approximate the dynamics of gas of particles 
with attraction. The principle difference is that the dynamics of brownian particles 
(RKSE) is diffusive (originates the overdamped motion with random force), while the 
dynamics of Bose atoms is the quantum analog of Newtonian mechanics. In both cases 
collapse occurs if the number of particles is large enough to cause attraction overcoming 
either quantum pressure (NLSE) or diffusion (RKSE). 

3. Self-similar collapsing solution of the 2D reduced Keller-Segel equation 

2D RKSE (0Q) is invariant under the scaling transformations p(r,t) — > -^p^r, j^t), 
c(r, t) — > c(j-r, -ht) for any L(t) = L = const > 0. Similar property holds for NLSE. 




(21) 



N c = 8tt 



(22) 
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The 2D RKSE has a static, radially-symmetric solution 

8 

(23) 



2\ 



c = -21n(l + r 

which corresponds to the critical number of bacteria, N(p ) = N c = 8n. This property is 
another striking similarity with the ground state soliton solution ijj = R(r)e lt , R(r) > 
of NLSE containing exactly the critical number of particles, P c = J R 2 dr. 

Assume that collapse is centered at r = 0. Then the solution of RKSE in the limit 
t — > t c approaches a radially-symmetric, self-similar solution. The self-similar solution 
has the form of the rescaled stationary solution (T23]) with a time-dependent scale (the 
collapse width) L{t): 



L(t) 2rU \L(t) 



c(r, t) = c 



r 



Lit) 



(24) 



The scale L(t) approaches zero for t — > t c . 

To describe the radially-symmetric solution we introduce the new dependent 
variable m as follows, 

m(r,t) = -!- / p(r',t) dr', (25) 



271 



r' <r 



which allows us to rewrite RKSE as the closed equation for m [4]: 

d t m = rd r r~ 1 d r m + r~ 1 md r m. (26) 

Here, m(r,t) has the meaning of the mass (the number of bacteria) inside the circle of 
radius r (up to a factor 2tt). Boundary condition for m at r — > oo is simply related to 
the total number of bacteria: m\ r=OQ = N/(2tt). In contrast to RKSE, Eq. ( 1261) is PDE 
for m. This simplification is possible only for radially-symmetric solutions of RKSE. 
In terms of m, the steady state solution ( 1231) of RKSE takes the following form: 

4r 2 

m » = (27) 
and the self-similar solution ( )24l) becomes 



itself simile 



i + y 2 



y = T (28) 

The boundary condition at infinity gives the critical number of bacteria, 
2nm se if similar | ^ —7-8^ = const, It also indicates that bacterial collapse is strong 
because the number of bacteria trapped within the collapse is nearly constant. 
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Assuming a power law dependence L(t) oc (to — t) 13 of the collapse width in the 
self-similar solution ([28]) one concludes that all terms in Eq. fl26|) are of the same 
order provided (3 = 1/2, which is similar to NLSE where also the collapsing width 
oc (to — t) 1 / 2 . Like for NLSE, the self-similar solution ( )28|) is not an exact solution 
of Eq. ( 1261) . To account for the difference, it is necessary to consider the logarithmic 
correction to L(t) ~ (t - t) 1/2 : L = (t - t) 1/2 f(\n (t - £)), where /(In (t - t)) is a 
slow function compared with (i — t) l l 2 . This slow function comes from the nearly exact 
balance between linear and nonlinear terms of RKSE (between diffusion and attraction). 
The same slow function allows to introduce a small parameter a, defined in (jSJ), which is 
a slow function of (to—t) 1 / 2 compared with L. The balance between linear and nonlinear 
terms of RKSE improves with decrease of a — > 0. 

Based on the analogy with the critical NLSE, we introduce in Eq. ( l26l) the new 
independent "blow up" variables [T6] : 



r 

y 



dt> ( 29 ) 



T 



These new variables transform Eq. ( 126]) into the equation for a new unknown function 

(p(y,r) =m(r,t) (30) 

into the following equation 

d T ip = yd y {y~ l dy(p) + y~ x <pdy(p - ayd y tp, (31) 

where a is given by ([51) . The advantage of working in blow up variables is that the 
collapse occurs at r = oo instead of t = t c , so that the collapse time t c is eliminated 
from consideration. Also, the function tp has bounded derivatives. 

Figures EH,b shows that as t — > t c , the density p(r) grows near r = while the 
tail of p(r) is practically frozen for r > 3 (on a timescale of collapse). In contrast, the 
solution in the blow up variables is steady at y < 1 and is well-approximated by (1231) 
and (12^1) . as shown in Figure [3b, d. It is also seen that the deviation of solution from 
(f23[) moves away from the origin y — as t — > t c . 

Based on our assumption that a is a slow function, it is natural to look at the 
solutions of Eq. (j3T!) in the adiabatic approximation where one can neglect r-derivative 
in the left-hand side (lhs) of Eq. [3TJ Then, assuming that |a| 1, one can expand the 
solution of ( l3~TT) in powers of a starting from (1281) for the power zero. Unfortunately, 
the term —ayd y ip grows with y and violates the expansion for large y. So, the adiabatic 
approximation can only work locally and is restricted to not very large y, a situation 
which is familiar from the analysis of collapse in NLSE. This however does not create 
a problem because the behavior at large y does not affect the self-similar solution near 
zero. 




Figure 3. The spatial dependence of the density p (panels a,c) and the mass m 
(panels b,d) at different moments of time for the simulation with N = 1.0250iV c . In 
the top row (panels a,b), the data is shown in simulation coordinates. In the bottom 
row (panels c,d), rescaled density, L 2 p(y), and mass, <p(y,r) = m(r,t), are shown as 
functions of rescaled radius, y = r/L. In panel (a), notice the growth of p(r) near 
the origin and a nearly steady tail. In panel (c), notice the convergence to the static 
solution ([2B")) . in the growing neighborhood of y = 0. In loglog scale the deviation 
from that static solution has the form of a bump. The bump moves away from the 
origin as t — > t c . 
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It is convenient to present a general solution of Eq. (13 ip in the following form 



A 2 2 

r) = + et^ 2 -^—v(y : r), (32) 



where v(y, r) includes all corrections with respect to the self-similar solution fl28|) . Here, 
the factor e^ y2 (which plays a role of a gauge transform) is inspired by a somewhat 
similar factor e~ l ^ y2 in the self-similar solution of NLSE [31]. However, the absence of 
—i in the exponent makes RKSE case quite distinct from NLSE case. 
Substitution of ( 13"2"j) into ( 13~Tj) gives the following equation: 

d T v + l a v = F. (33) 

Here 

X 1 ^ 3 „ 8 rfl n 2(3 -| . . 

C = --d y y 3 d y - + -y 2 - 2a + — — , 34 

?T (l+?r) z L 4 l + y^ J 

is the linear operator corresponding to the linearization of ( 131 j) with respect to (I28p . 
The right-hand side, 

F = -^V, - ^e"- 2 /4 + V/4 + 2 ^ + ^ e ^/4 5 (35) 



y 2 + l 2(y 2 + l) (y 2 + l) 2 y 2 + l 

is responsible for all other terms. These other terms include terms nonlinear in v, 
inhomogeneous terms, and linear terms. Generally, F cannot be zero because (|28|) is 
not an exact solution of ( l3~TT) for nonzero a. Notice that up to now we have not made 
any approximations, so Eqs. ()32l)-(l35l) are equivalent to Eq. ( l3~Ti) . 

The advantage of the definition (132]) is that the operator C a = —^d y y 3 d y + V(y) 
has the form of the radially symmetric Schrodinger operator in spatial dimension four 
(D = 4) with the potential 

nv) = -77^ + & - 2a + tt^] • ( 36 ) 

(1 + 2/ ) 4 l + y 1 -* 

It means that C a is the self-adjoint operator with the scalar product 

oo 

<<M>= / ^{y)cj>{y)y^dy. (37) 



The potential V(y) — > oo for y — >• oo, which ensures that C a has only discrete spectrum. 
This allows us to expand arbitrary v in a discrete set of eigenfunctions of C a : 

v = ci^i + c 2 ^ 2 + c 3 ip 3 + (38) 

where c\(r), C2(r), ... are r— dependent coefficients of the expansion (below we often 
omit argument r for brevity), ipjii/) are the eigenfunctions of £ a , 

Caljjj = Xjlpj, (39) 
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and Xj are the respective eigenvalues. The eigenvalues are ordered starting from the 
lowest eigenvalue as Ai < A2 < .... All eigenvalues are real and non-degenerate as 
discussed in the next section. 

Note that the use of the scalar product (|37|) (which corresponds to the radially 
symmetric Schrodinger operator in D = 4) is simply an auxiliary mathematical trick, 
which is effective because the operator C a is self-adjoint with this scalar product. We 
remind that all solutions obtained below correspond to RKSE ([I]) with D = 2. 



4. Spectrum of linearization operator 



The eigenvalues ( 1391) of the linearization operator C a are given by the following implicit 
expression, 



A + 2a 
2a 



In ^ 

a 



A_ 

2a 



+ K 



1 + a 1 / 2 In 



1 



(40) 



as it was proven in Ref. [15] using a rigorous version of the method of matched 
asymptotics. Here K := In 2 — 1 — 27, while ^ is the digamma function, defined as 
\l/(s) = 4- lnT(s), where T(s) is the gamma function. 

Solving (14"0]) for A gives the spectrum of C a , starting from the lowest eigenvalues, 
as follows 



Ai = a 



-2 + r^ + 2(l+ 7 -ln2 N 



A 2 = a 



In 



2 , , 1 

2(2 + 7-ln2 



(ln± 



l\2 



2{K + 1 f 



7T 



(ln± 

v a 



l\3 



(ln± 



L\2 



-4 + 2(^ + 7)^ 



7T 



o 

2 

o 



OO 4 
1 



(lni 



1^3 



0^5 



1^4 



(41) 



A, = a 2 



^ T + (5 + 2 7 -21n2) 7U \- 



lni 



+ 



2(1 -3A' -37) + 2(K + rf 



21 



(lni 



l\3 



o 



(lni 



1^4 



Eigenfunctions can be also approximated from the method of matched asymptotics. 

In Section [5] we need to calculate multiple integrals which involve ipj. For this 
purpose, it is more convenient to use the variational approximation for eigenfunctions 
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obtained in Ref. 



V'3 



-ay 2 /A 



1 1 

ay ay 



In (1 + y 2 ) ) e~ ay2/4 



ay 



tt 2 \xi\ -12 + 7T 2 (2 + 7 -ln2) 



7T 



12 



7T Z 



12 



(42) 



-ay 2 \n{l+y 2 ) 



12 



7T 2 - 12 



4 



In - -3-7 + In 2 
a 



24 



7T 



12 



-aj/ 2 /4 



where ipj means the variational approximation to ipj, j = 1,2,.... We estimate the 
accuracy of the variational approximation by calculating the variational approximation 
for the three lowest eigenvalues Ai, A 2 , and A 3 as A,- = Shj^d^l, j = 1,2,3. These 



scalar products involve the calculation of integrals of the type described in Appendix B 



Expansion of the resulting expressions for integrals Ai, A 2 , and A3 in inverse powers of 
In 1 agrees with exact results f )4TT) up to order O ^ ln i) 2 ) f° r ^2 and up to order 

O f^r) f° r ^3- This is the best result we are able to achieve with the variational 
approximation. This accuracy will be however sufficient to obtain ([7]). 



5. Amplitude equations 

Similar to fl38|) . we expand v from (133]) in a set of approximate variational eigenf unctions 
t/jj, j = 1,2, ... as follows, 

00 

v = ^2 c $v ( 43 ) 

where Cj(t) are the coefficients of the expansion. In this Section we derive a set of 
amplitude equations for Ci(r), c 2 (r), . . . from fT43|) which provide a solution of Eq. fl33|) . 
We solve the amplitude equations exploiting the fact that, at the leading order in a, the 
solution of ( 13"T|) is given by f ]2"g]) . (We used that fact in the definition of fl3"2"]) ). We expand 
all expressions below in integer powers of the small parameters a and ^p-, keeping the 

a 

lowest nontrivial order of a and several orders of r-nr- 

In - 

a 

We assume the approximate orthogonality of the variational functions, 

(M 3 ) = 0(a)U i \\\\^ j \\ for i^j, (44) 

where := (ipi^i) 1 ^ 2 ^ i — 1,2, . . .. Then, the scalar multiplication of (133]) onto ipj 
(with the scalar product (|3Tj) ) results in 

(tpj,d T v) + (i!j,C a v) - (tpj,F(v)) 

00 00 (45) 

i=i i=i 
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Here, we neglect corrections from nonexact orthogonality (144 j) because, as we show later, 
these corrections are of the next order in a when compared with other terms in fj45|) . 
In this section, all calculations of scalar products for ( 145]) are based on integrals defined 
in Appendix B for the variational functions ( 142]) . For instance, the direct calculation 



for the variational functions ( 1421) gives the following expressions: 

||^i|| 2 = -321na + 32(-l - 7 + ln2) + 0(alna) , 

||^2 f = 32(lna) 2 + 32(1 + 2 7 - 2 In 2) In a, 
16 

+ — [7r 2 + 6(hi2- l)ln2 + 67(7+ 1 -2 In 2)] + 0(alna), 
o 

I, 7 ,12 aA n , 2 , 32 (-108 - 487 + 137r 2 + 4 7 7r 2 + 481n2-47r 2 ln2) l 

H^ll =64(lno) + (-12 + lna (46) 

+ ^32 h 7 2 (-12 + vr 2 ) 2 + vr 4 [23 + In 2(-13 + 2 In 2)1) 

(-12 + 7T 2 ) 2 L v ' 

- 247r 2 [13 + ln2(-ll + 2 In 2)] + 144(-1 + In 2)(-7 + 2 In 2) 
-7(-12 + tt 2 ) [108 -481n2 + 7T 2 (-13 + 41n2)] + 0(a(lna) 2 ). 

We assume (based, e.g., on numerical simulations in [HI [8] and following Ref. [14"] ) 
that a is the adiabatically slow function of r: d T a <C a 2 . As mentioned above, we 
expand all quantities in the small parameters a and -^r (it is also seen in Appendix B 



that all integrals involved in (145|) expand into these parameters) keeping only a 
leading order in a and many enough terms in powers of ^V. Then the adiabatic 

a 

assumption d T a <C a 2 requires d T a = a 2 (j^tj • We introduce a normalized function 
«t : = ^j^t = 0(1) + 0(a). This allows to write d T a as an expansion in inverse powers 

a 

of In i only: 

d T a = a 2 -^ T ~a T , ~a T = 4 0) + S?> r^j + a?-\-^ + O (jr\^) , (47) 

where the coefficients hr and are 0(1) and do not depend on r in the adiabatic 
approximation. Note that the subscript r in these coefficients is not a partial derivative 
but rather indication that these are the expansion coefficients for a T . 

Assume that the expansion coefficients c\, C2, cs, ■ ■ ■ in f|43|) are initially 0(1). 
A series expansion of equations ( 145]) over small a, using Eq. ( 1471) and dividing each 
jth equation by ||^j|| 2 , together with (|46j) . result at the leading order in the following 
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expressions 



d T c\ + a — 2ac\ + O ( j = 0, 
d T c 2 + O (°) = 0, 
d T c 3 + 2ac 3 + O (^j) =0, 
d T c 4 + 4ac 4 + ( p\ J = 0. 



(48) 



Here, the terms 2a(j — 2)cj, j = 1, 2, 3 originate from eigenvalues for ipj (see Eq. f )4T]) ). 
while the term a in the first equation comes from the scalar product of ipi with the second 
term in the right-hand side of (|35|) . Also the contribution from d T ipj = (d T a)d a ipj, j = 
1,2, . . . is included into 0(. . .) term. It follows from Eqs. (148j) that the coefficient C3 
initially decays exponentially (because a > 0) until it reaches the adiabatic, quasi-steady 
state with C3 = O (j^t)- O ur conjecture is that the other coefficients, 04,05, . . ., also 
decay exponentially (they correspond to the larger values Xj, so that they are assumed 
to decay as Cj oc exp [— 2a(j — 2)r], according to the linear terms in (133|) ). The lack 
of explicit expressions for ipj, j > 4 does not allow us to prove this statement. We 
conclude that, after an initial transient, the coefficients C3, C4, . . . reach the adiabatic 
state with their values 



a 



c 3 , c 4 , . . . = O ( — r ) . (49) 



Below we assume this adiabatic state. 

In the first equation of ( 1481) we assume that 



c i = S + ( T^T ) 



2 \ ln± 

v a 

to avoid exponential growth of ci in r. (Such artificial exponential growth would result 
in error in estimating t c .) 

We have now a freedom in selecting c 2 , and we choose it so that v — > for any |/ 
as a 0. According to (1421) . ^(y)!^ = ^2 (2/) | y=o = 8 so we set 

In this case, ciipi + c 2 ^ 2 = 0(a) for y = 0(1), i.e. w in fl32l vanishes with a — > 0, as we 
expect from the self-similar solution ( 128]) . 

Equations. (149|) . ( )50|) . and (|5T]) justify the adiabatic approximation, which means 
that the coefficients ci, c 2 , C3, C4, . . . depend on r only through a, and one can expand 
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Cl 



c 2 



c 3 



fe=i 

oc 



1 - . 1 



fc=l 

oo 



(52) 



k=l 



where the expansion coefficients 4 = 0(1) for any i, j; the coefficients do not explicitly 
depend on r in the adiabatic approximation. 
It follows from ([52} and flS} that 



d T Cj —Old. 



r lni/ ~V(ln±) 3 

Similar to derivation of Eqs. ( 148|) . we now perform a series expansion of equations 
into small a (but in contrast to the derivation of Eqs. (I48p we proceed to the higher 
orders of expansion) using Eqs. (j4"6"|) . ( l4"T|) . f )52|) to obtain following equations: 



O 



J = 1,2,3, 



(53) 



d T ci + 



+ 



ln± 

a 



(i) 



(ln± 



1^2 



2eq 



+ 24 1) -24 2) + 24 1) -4 0) 4 1) + 24 1) 



o 



l\3 



(54) 



<9 T c 2 + 



ln± 

a 



-1 - 



+ 



(ln± 

V a 

a 



l\2 



(ln± 



1^3 



-(2) 



24 1} + ^(d^ - 24 1} ) - 7 + In 2) 



1 _ ^_ + srfW + (2 + 4°))4 2) + 64 1} + 4S(°)4 1) - 2fi«4 1) 



2al 0) 4 2) 



3 + 



7T 2 . 244 0) 4 1) 



6 -12 + 7T 2 

21n2 + 7 (-2-7 + 21n2) 



(In If + 4 ij (4 + 4 1J + 2 7 - 2 In 2) 
a 



O 



0. 



(55) 



d T c 3 + 



ln± L 

a 



2d 



(i) 



+ 



(ln± 



l\2 



-l + 2(l+4°))4 1) + 2rf 
+ 



1>3 



(56) 



Logarithmic scaling of the collapse in the critical Keller-Segel equation 



20 



Here we have neglected the expansion coefficients Cj for j > 3 by setting 04 = 05 = 
Cq — . . . — 0. In Equations fl54|) - fl56|) we keep the necessary number of orders in j-^- to 

a 

obtain the closed expressions for the expansion terms in (1471) . Equations (l54j) - (!56l) can 
be viewed as the compatibility conditions which ensure that expansions (J52l) and (|47j) 
are correct, so that a is indeed the adiabatically slow variable. 
It follows immediately from Eq. ( 156]) in the order that 

a 

4 1} = 0, (57) 
and from Eq. fl55|) in the order that 

a 

a? } = "2. (58) 
Then, from Eq. ( 154)) in the order we obtain 

a 

d? = -\ (59) 
Using Eqs. ( )55|) and (I57j) -( l59|) we obtain in the order ^i-^ that 

4 1} = -2 - 2 7 + 2 In 2. (60) 
Using Eqs. ( 154|) and ( 1ST)) -( 160]) we obtain in the order ^pr^ that 

4 2) = ^(-2 + 44 X) -7 + ln2). (61) 
Similar, using Eqs. ( 156]) and ( 15T|) -( 16T?]) we obtain in the order p^r^ that 

4 2) = 5. (62) 
Equation (|55|) in order , 3 requires also to take into account d T c 2 which is given by 

a — 44f +0 (4f)' (63) 

according to (1331) and ( 1521) . 



Using Eq. ([55} in order and (J57D-([60D, (P|) . (pj) we obtain the closed 

expression 

4 2) = 2(ln2) 2 + 41n2 + 7(-4-2 7 + 41n2). (64) 



Here, the unknown coefficient has been cancelled out identically. 

Equations (H7J, flSED, dHOD, and flBl} result in closed ODE (J7J) for a, which is 
the first main result of this paper. Figure H] shows d T a as a function of a for RKSE 
simulations with different initial conditions (the same initial conditions as in Figure 
[1]). Notice that after an initial transient all curves collapse to the single curve given 
by Eq. ([7]). This suggests that we can use the proximity of numerical curves to 
the analytical curve as the criterion for selecting to i n equation ( TlOl) . In Figure HJ 
we used the values of to defined for each initial condition as the time t = to when 
the relative difference between numerical and analytical curves reduces down to 20%. 
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a 



Figure 4. Dependence d T a{a) extracted from RKSE simulations shown in Figure Q] 
(thick solid, dashed, short-dashed and dashed-dotted lines, respectively). The thin 
dotted line represents the analytical dependence d T a(a) from Eq. (0, with neglected 
0{. . .) term. Notice that after the initial transient all numerical curves collapse on 
the analytical curve. The arrows point to the locations where the relative difference 
between the analytical and numerical curveds reduces to 20%, the criterion for selecting 
t in Eq. (flQ|) and in Figure [TJ 

Arrows in Figure H] point to locations (a(to) , d T a(to)) satisfying this criterion. For the 
simulations of Figure [1] we obtained t = 7.2125 . . . , t = 4.5879 . . . , t — 3.3257 . . . , 
to = 2.5528... for N/N c = 1.0250, 1.0375, 1.0500, 1.0625, respectively. Also in these 
cases t c = 8.12305 . . . , t c = 5.32533 . . . , t c = 3.94247 . . . , t c = 3.12039 . . . , respectively. 

The dashed-dotted curves in Figure [1] are only weakly sensitive to the choice of 
to < t c , provided to is chosen later than the time specified by the 20%-difference criterion. 
For instance, if we choose t based on 10%-difference criterion (instead of 20%), the L{t) 
curves in Figure [U would change by < 5% which is within the relative error of these 
curves in comparison with the numerics (solid curves in Figure [T|). 

It also follows from Eqs. ( J52l) . (JBTl) . ( |59l) . and (J62l) that the expansion coefficients 
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in f )43p are given by the following expressions 

— H +0 (4?)- (65) 

Thus, the coefficient C3, which corresponds to positive eigenvalue, is of a lower order 
compare with c\ and c%. We expect the similar to be true for all coefficients C3, C4, . . .. 
Also, the coefficient is undetermined in our approximation order. We expect 
that it might depend on initial conditions. We conclude that the self-similar solution 
(12 8p is stable with respect to radially-symmetric perturbations, and the leading order 
corrections to it are determined by a linear combination v ~ c\ipi + C21P2, where c± and 
c 2 are given by (|65|) . 



6. Blow-up rate of self-similar solution 

In this section we solve ODE ([7]) together with ([8]) and (Q to derive the blow-up rate ()3]). 
Integration of Eq. (0) from an initial value To to r gives 
1 



a 



ln I + M+ ^ 

a in- \[\n-) 2 



2(r-r ), (66) 



-TQ 



where M = f - 1 = -2 - 7 + In 2 and b = | + Mi = 1 + nl as i n if we look at 

Eq. flB^l) as the implicit expression to determine a(r) then it turns into a remote relative 
of the Lambert W-function. Such implicit expression can be solved for a assuming r 3> 1 
by iterations as follows: 

1 L 2 L 2 2 L 2 M —lb + 2M + M 2 — 2ML 2 _ . , 



where L x := In [2(r - r*)], L 2 := In In [2(r - r*)], and 

a = a(T ). (68) 



r* = r - — ( In - + M 



CC=CCQ 



2a \ a In - 

At this point, one can proceed in qualitatively similar way to Ref. [H] to determine 
L(t). However that way of calculation results in a slow convergence of the asymptotic 
series for L(t) with the increase of r. We choose a different path. Our goal is 
to start with Eq. (166 p . to carry as many steps of exact transformations as possible, 
and to perform asymptotic expansions as late as possible. Here and below we abuse 
notation and use the same notations for all functions with the same physical meaning, 
independently of their arguments: L = L{t) = L{r) = L(a), r = r(£) = r(L) = r(a) 
and a = a(t) = a(r) = a(L). Similar, for initial values L = L(t ) = L(t ) = L(ao), 
tq = r(t ) = t(L ) = r(a ) and a = a(t ) = a(r ) = a(L ). 
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(69) 



We integrate (169]) in r between To and r, using the integration by parts, to obtain 

l r T r a 

— In — = / a(r) dr = ar(a) — a r(a ) — r da 

L$ J Tn Jan 



r - t J a - r - r a 



(Y — r*) da. 



(70) 



To evaluate integral over a in (170]) explicitly we use r(a) from (166]) with (168]) and 
obtain 



In — = - 

L 4 



n - ] n — 

a 



M + l 



+- ( In In In In — 

a ao 



2 \ In i In 

a ao 



In In — 

a a 



1b i 



(71) 



Note also that the term O ^^rj m dZl]) originates from the next order term O (^j^t 

in Eq. ( 166]) . Formally, in Eq. flTT]) . the terms O ( ) and are of the same order. 
Yet, our numerical simulations indicate that ^p-x term improves accuracy of the analytic 



21n- u 

approximation, so we keep this term in its explicit form. 
We introduce new variables, 

I := In — and / :=ln— , 
L Lq 



as well as define 



r = / _-(ln-) - 

ao 



M+l 1 



In I In In — 

a 2 



which allows to rewrite (ITT]) as follows: 

2 



Z-Z* 



In 



M + l 1 b 
In - + - 

2 a 2 



In In 



a In - 



a In — 

ao 



O 



i J ' 



(72) 



(73) 



ln± 



(74) 



We now solve Eq. (1741) for ln^. Instead of doing straightforward iterations, we 
neglect the terms |(. . .), O ^r^rj in Eq. ( 174]) and solve the remaining part of the 
equation, F 2 — 2Y — V = 0, exactly: 



Y = 1 + VT+V, 



where we define 



Y :-- 



lni 

a 

M + l 



and V :- 



(M+l) 



(75) 
(76) 
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with Yq being the leading order approximation to Y, such that 

Y = Y + 6Y. (77) 

To find 5Y as a function of V, we represent SY through the formal series 
SY = J2™=i & ~Y^ ( w ith Yq given by (175]) ). We use this series together with (T75|) -(f77I) to 
perform a series expansion of Eq. (1741) in inverse powers of Yq. It allows to determine 
the coefficients 5Y n recursively at integer inverse powers of Yq starting with the power 

bu\-(l+M)Y ] 



zero. In particular, the zero power gives YLi 



(l+M) 2 



Note that the double 



logarithm In In - in ( I74p also needs to be expanded. All together it results in 



b In 



5Y 



[1 + M)Y -6(1 + M)ln -{1 + M)Y 



[i + m) 2 y 



1 + MfY 2 



+0 



Y 2 



-0 



(InFo) 



Y 3 



Here, similar to (]71j) . we keep the term — 



(78) 



3-^, even though this term is of the same 



order as O term. Here, M + 1 = -0.884068... according to (TTUl) . Note, that 

instead of performing an expansion in inverse powers of Yq, one can simply do it in 
inverse powers of V 1 ! 2 . This, however, would result in a slower convergence for moderate 
(V > 1) values of V. 

We rewrite OH]) as — = dt, and integrate it between time t c and t: 



dt' = t r -t 



-^—AJJ, 

a(L') ' 



(79) 



where following 
(M 



exp 



I* + 



2) 

l) 2 



and (T76]) we can represent L through V as L 



-V 



The dependence a(L) in ( 175^) follows from ( 172^) -( 175|) . 



Switching from integration over L to the integration over Yq in ( |79l) we obtain: 



t, - 1 



exp —2 



r + 



(M 



x exp -(1 + M)F ' 



-6(1 + M) In 



6 In 



-[(^0- 
M)Y Q > 



1] 



(M + l) 



;i + m)f ' 



;i + M) 2 r ' 5 




(80) 



Here, the integration cannot be carried explicitly. Instead, we use the Laplace method 
(see e.g. [Hil HZ]) to evaluate the integral asymptotically in the limit Yq ^> 1. We 
introduce in Eq. (|T0|) a new integration variable, 



2 :— F — Yq, 
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(M+ l) 2 
t c — t = exp 



-21* - ( M + 1 )V 2 + (M + 1)MF 



e^'^dz, (82) 



where 



S{Y ,z) = -{M + l) 2 z 



(M + iy 



z 2 + (M + l)Mz + ln(y + z - 1) 



6 In 



+ 



-(l + M)(y + z) -6(1 + M)ln _(l + M)(y + z) 



+ 6 



y (i + M)(y + ^) 



y (i + M) 2 (y + z) 2 
+ o 



y (y + ^) 2 



To use the Laplace method for asymptotic expansion of the integral in 
start with the following general expression, jlOJ H7] : 



/ e YoS(z,Y ) dz _ e Y S(0,Y ) Cn Y~ n - 1 
n "=0 



(83) 



we 



34) 



with 



■1) 



n+l 



5 



z=0 



S'{z,Y ) := -^S(z,Y ) 



55) 



Taking into account two leading terms in ( 1541) . we obtain from ( 1521) . ( 1531) . ( 1541) . and ([HS 
the following expression: 



(M+ l) 2 
t c — t = exp 



-21* 



(M+iy 



y 2 + (M + 1)MY + ln(y - 1) 



x exp 



b In 



(1 + M)y -6(1 + M)ln -(1 + M)y 



1 



(i + M)y 



i + 



(i + M) 2 y 2 



+ 6 / J 



M 



+ 



M 2 



+ 



/i n y 



(M + i)2y (i + M)y (i + M) 2 y 2 V ^o 3 

We now define a large parameter 

x := ^/-21n/3(t c - t) - M, 

where 



(86) 



57) 



2 exp < 21* - — }■ 
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We multiply both the lhs and the rhs of (]55"]) by f3 from 0881) and take logarithm 

2 

from both sides to obtain -|- on the lhs. We solve the resulting equation for Yq by 
assuming the asymptotic form, 



Y = b^x + 



n=0 



X r ' 



(89) 



and performing a series expansion of both rhs and lhs of that resulting equation in 
inverse powers of x. The coefficients 6_i, 6 3 are determined recursively giving 

l-ftlnsc -\ -2M + b(\nx + 2M\nx- 1 



Yn 



M + l 



X"' 



O 



(In a;) 



(90) 



x° J \ x* 

Note that the choice of the factor exp j-^j in fl55J is somewhat arbitrary (the lhs 
and the rhs of ( 1561) can be multiplied by an arbitrary positive constant). The factor 
exp | — is chosen to speed up convergence of ( 1901) for Y > 1, i.e. for L(t) < 1. 



Using (J72D, fllSD, and 0761) we obtain 



L(t) = exp 



{M + l) 
4 



<Y 2 -2Y ) 



(91) 



Equations (JTJ), (JTJ) ([57]), flSSJ, O, and ([21]) give the closed expression for as 
a function of t c — t and the initial values L = L(to), a$ = —LL t \ t=to . To make 
the comparison with the old scaling ([6]) more transparent, we plug the expression 
for Yq from ([90]) into ([9T]) and perform a series expansion of the expression in the 
exponent into inverse powers of x, obtaining the final expression ( [TO]) . Note that the 
first term in the exponent of the first equation in ([TO]) can be rewritten through x as 

~ V ~ — — X+ 2 I ■ Thus, Eq. ([TO]) includes terms of orders x and x°. 

The error terms O (^-) in ([9T?]) and O (^j in ( 1101) result from the error term 
O y ^rji j m ^' c l- ©• however, chose to write down explicitly the terms of the 

same orders, oc in ([90]) and oc \ in ([TO]) . These terms are independent from the error 



term O 

m am. 



(ln± 

v a 



of Eq. OTJ) . Next order error terms are O 



in (jSnP and O 



(lna:) ; 



7. Numerical simulations of RKSE 



In our numerical simulation we evolve Eq. ([26]) . written in terms of the mass of bacteria 
m(r,t) within the circle of radius r as defined in ( [25]) . The density, p(r,t) = £^r, and 
other quantities characterizing the evolution of the collapse are computed from the mass. 
To find the width of the collapse, we assume that the solution has reached its self-similar 
form given by Eq. ([241) . Then, the collapse width can be estimated from the density at 
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m 




Figure 5. Schematic representation of the discretized solution and the grid structure. 
Three subgrids are shown. The subgrids closer to the center of the collapse have finer 
resolution. The data at black points are evolved by the discretized Eq. (f26|) . the data at 
white points is copied from neighboring subgrids (the copying is shown by arrows), the 
data in gray points is interpolated from neighboring points using 6 order polynomial. 



the center as L = (Ip^o)" 1 ^ 2 - To compute the slow parameter, a, we differentiate L(t), 
as in OH]). The self-similar time, r, is found by integrating L(t) according to Eq. Q. 

A typical solution for m(r, t) is shown in Figure [3b. The spatial extent of the 
collapse is marked by the large gradient of the solution near the center, which becomes 
even larger and moves even closer to the center as time progresses. This requires special 
treatment to ensure that the solution remains well-resolved. 

The results presented in this paper are obtained using an adaptive mesh 
refinement (AMR) technique [HJ [31], complemented with the fourth-order Runge- 
Kutta time advancement method. Our spatial domain, r e [0,r max ], is divided into 
several subdomains (subgrids) with different spatial resolution. The spacing between 
computational points is constant for each subgrid, and differs by a factor of two between 
adjacent subgrids. The rightmost subgrid, farthest from the collapse, has the coarsest 
resolution; the spatial step decreases in the inward direction. 

The grid structure adapts during the evolution of the collapse to keep the solution 
well resolved. When a refinement condition is met, the leftmost subgrid is divided 
in two equal subgrids. Then, the new leftmost subgrid is refined; that is, additional 
computational points are placed halfway between the existing points. The values at 
the new points are obtained with sixth-order interpolation. The condition for grid 
refinement comes naturally from properties of the self-similar profile. The density at 
the origin increases by a factor of 4 as the width of profile decreases by a factor of 2 
according to Eq. ([3]). To keep the effective number of grid points per L within desired 
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limits, we use the increase of the maximum density by factor of 4 as the condition for 
refinement. 

In the interior of each subgrid, the spatial derivatives are computed using fourh- 
order central differences on the five point stencil. At the subgrid boundaries, the data 
are copied between subgrids to fill in values at "ghost points", as shown by arrows 
in Figure [7J Notice that communication between subgrids is going in both directions: 
in AMR terminology, data from the fine subgrid is restricted to the coarse grid ghost 
points, and data from the coarse grid is prolongated to the fine grid ghost points. The 
data between points of the coarser subgrids needed for finer subgrid ghost points are 
obtained by sixth order interpolation. The left ghost points of the leftmost subgrid are 
filled using reflective boundary conditions. The point r = is treated in a special way 
because of the singularity in the rhs of Eq. ( |26|) . Expanding m(r, t) in a power series in r 
at the origin and using the definition (1251) we obtain that m(r, t) = p ^°^ r + 0(r 4 ). This 
is also consistent with the series expansion in r of rhs of Eq. (I26p . Thus in the spatial 
discretization we set m(r — 0, t) — 0. The right ghost points of the rightmost subgrid 
are filled with the data from the last point. We found the right boundary conditions 
to be very forgiving, which is not surprising considering that the mass approaches a 
constant, as r~ 3 , when r — > oo. 

The solution on all subgrids is evolved with the same timestep, At = CcFLh 2 , 
where h is the spatial step of the finest grid and Cqfl is the constant. We typically 
used Ccfl = 0.4 but also tested a convergence for smaller values of Cqfl- 

We use two kinds of initial conditions. First kind is the Gaussian, 



which implies 

a 1 

where a and A are the parameters of the initial condition. Second kind is the modified 
stationary solution, m\t=o = Amo(r), where mo(r) is given by Eq. (1271) . Both types of 
initial data result in similar dynamics for the same values of N = lixA. The simulations 
presented in this paper were performed for Gaussian initial conditions, with a — 1 and 
A = 4.1, 4.15, 4.2, and 4.25. The initial grid was comprised of ten subgrids; the finest 
subgrid had 400 points while all other subgrids had 200 points each. The size of domain 
was set to r max = 1600L . 

We have verified the AMR code against an independently developed, uniform grid 
code with an adaptive spatial resolution and an adaptive time step. Similar to the 
AMR code, the uniform-grid code evolved Eq. (liZo]) using fourth order Runger-Kutta 
integration in time. The spatial derivatives were computed spectrally using the FFTW-3 
library [49]. Since Fourier transforms require periodic boundary conditions, the spatial 
domain was extended tore [— r max , r max ] with sufficiently large r max (about 50L ). As 
in the AMR code, we set the value of the mass to zero at the origin point to avoid the 
singularity in the rhs of the equation. Although uniform in space, the grid resolution 
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was refined at times when the maximum density increases by a factor of four. The new 
grid had twice as many points with values computed by spectral interpolation. We run 
the uniform-grid code at Cqfl = 0.2. 

Although the uniform-grid code was useful for cross-comparison, it was significantly 
less efficient than the AMR code. Typically, we run the uniform grid code until the peak 
density reached ~ 10 5 , which required 32, 768 gridpoints with grid resolution ^ < h < j. 
On the other hand, in the AMR simulations presented here, the density reached ~ 10 17 
on approximately 12, 000 total gridpoints, with ^ < h < ^ resolution on the finest 
subgrid. 

8. Conclusion and Discussion 

In conclusion, we studied the collapsing solution of the 2D RKSE, Eq. (CQ). To leading 
order, the collapsing solution has the self-similar form ((3]), characterized by the scaling 
L(t). Our analysis of the dynamics of perturbations about the self-similar form allowed 
us to find the time dependence of the width of the collapsing solution given by the 
scaling f fTUj) . The analysis of the perturbations is performed by switching to independent 
"blow up" variables ( 1291) . and an unknown function (1301) . In the blow-up variables, the 
analysis of the dynamics of the collapse reduces to the analysis of the perturbation 
about the static solution (|28|) . The analysis exploits the slow evolution of the parameter 
a, defined in flEl), which originates from the leading order scaling L{t) oc (t c — t) l l 2 . 
After applying the gauge transform (1321) . we expanded the general perturbations in 
eigenf unctions of the self-adjoint linearization operator C a about the static solution 
and derived the system of amplitude equations. We solved these amplitude equations 
approximately to obtain ODE ((7|) for a(r). We solve Eq. (0) asymptotically in the limit 
t — > t c , together with (IE|) and Qj to obtain the scaling ( [TO]) . 

We found that both ODE ([7]) for a(r) and the scaling ([TO]) for L(t) are in excellent 
agreement with numerical simulations of RKSE. We compared the scaling ( [TO]) with the 
previously known scalings (jl])-© and showed that scaling ([5]) is the correct asymptotic 
limit. However, this limit dominates only for unrealistically small values L < 10~ 10000 . 
In contrast, the scaling ([TO]) agrees well with simulations for a quite moderate decrease 
of L(t) compared to the initial condition. E.g., Figure [2] shows that six-fold decrease 
of L compare with the initial value £(0) is enough to achieve the relative error < 7% 
between simulations and the scaling ([TO]) . 

We now discuss the limitations of the analysis of this paper. The analysis is 
exact until we derive the amplitude equation f[T5]) . At this point we must resort to 
approximation, because we can calculate only a finite number of terms in the amplitude 
equation. We approximate the eigenfunctions ipj of C a as i/jj, j — 1, 2, . . . through the 
variational analysis (see Section[4]). Such variational approximation itself does not create 
any obstacle because one can, at least in principle, expand the general perturbation 
about the static solution in functions ipj, provided they form a complete set on the 
space £ 2 ([0, oo), y 3 dy) which corresponds to the scalar product ( [37]) . However, the 
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variational construction of such functions turns out to be difficult for j > 3 (provided 
we aim to approximate the eigenvalues Xj of C a with high precision). This work is left 
for the future. In this paper we limit the analysis of the amplitude equations for ipj, 
j = 1, 2, 3 by setting Cj = for j > 3 in the amplitude equation. We found that already 
ip3 gives contribution only to the coefficient bo in ODE (JTj), i.e. to the highest order 
term . b ° 1)3 which we take into account. The other two lower order terms in rhs of ([7]), 

— j^t + n^i)2 1 are fuhy determined by ipi an d 4>2- We expect that taking into account 

the nonzero values of ipj and Cj for j > 3 might modify the value of bo in ODE for 
a(r) (and respectively modify b in (ITU!)). Such type of calculation presents a technical 
challenge and is left for the future. A potential route for such calculation might be the 
approximation of the eigenfunctions ipj using the matched asymptotic technique |15j . 
Using this technique, however, one faces the challenge of calculating the scalar products 
in the amplitude equation. It should be mentioned that the dramatic improvement of the 
accuracy of L(t) with the increase of the order of approximation (as seen in comparison 
of Figures [Tb,c with Figure [Ih) suggests that the essential part of b (and respectively 
b) is already captured in our scaling ([TO]) . 

Appendix A. Keller-Segel model of bacterial aggregation 

Bacteria and biological cells often communicate through chemotaxis, the process of 
secretion and detection of a substance called chemoattractant. Below we refer to 
bacteria and cell as synonyms. The chemoattractant secreted by bacteria diffuses 
through media. Other bacteria of the same kind detect it and move along its 
gradient. Thus the chemotaxis creates nonlocal attraction between bacteria. Bacteria 
are self-propelled and, without the chemotactic clue, the center of mass of each 
bacterium typically experiences a random walk. The motion of bacterial colonies is 
thus affected by the competition between random-walk-based diffusion and chemotaxis- 
based attraction. The macroscopically averaged motion of bacteria can be described 
by the Keller-Segel model (sometimes called the Patlak-Keller-Segel equation), see e.g., 
[H[2l[3l[ll[5l[6l[Tl[8l[9l [101 HU H21 [131 HH [HI US] and references therein: 

d tP = DV 2 p- Vpfe/oVc], (A.l) 
d t c = D c V 2 c + ap, (A.2) 

where p(r, t) is the bacterial density at spatial point r and time t, c(r, t) is the 
concentration of chemoattractant, D is the diffusion coefficient of bacteria (representing 
the random walk), D c is the diffusion coefficient of chemoattractant, a is the production 
rate of chemoattractant by bacteria, and the coefficient k > characterizes the strength 
of chemotaxis. 

The Keller-Segel model is a mean-field approximation of the behavior of a large 
number of bacteria, and can be derived from the dynamics of individual bacteria using 
macroscopic averaging over an ensemble of realizations of stochastic bacteria motion. A 
starting point of the derivation can be, e.g., the description of an ensemble of bacteria as 
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point-wise objects subject to a white noise force, as in Ref. [13]. Such description is most 
relevant to procaryotic bacteria like Escherichia coli which small rigid shapes. Another 
possible starting point is the description of the dynamics of eukaryotic organisms with 
randomly fluctuating shape, such as Dictyostelium amoeba [23 | |24" | [25]. 

If the initial density of bacteria is low, the bacterial diffusion typically dominates 
attraction and the density remains low. For instance, a typical time scale for the 
evolution of a low-density Escherichia coli distribution in a petri dish is about one 
day [5] (see Figure 3A in Ref. [5]). If the initial density is relatively high, attraction 
dominates, and bacteria aggregate (see Figure 3B in Ref. [5]). The typical time scale 
of such aggregation in experiments on Escherichia coli is several minutes [5]. Thus the 
aggregation has an explosive character compared to the evolution of bacteria outside the 
aggregation area. The aggregation is described by the "collapse of bacterial density" in 
the approximation of the Keller-Segel model flA.ip - f1A.2fl . 

The diffusion of chemoattractant is usually much faster than the diffusion of 
bacteria, i.e., D/D c <C 1. For instance, D/D c ~ 1/40 — 1/400 for the cellular slime 
mold Dictyostelium [50], and D/D c ~ 1/30 for microglia cells and neutrophils [5T| 152]. 
(Here, we refer to bacteria and cell as synonyms.) Thus Eq. ( 1A.2I) evolves on a much 
smaller time scale than Eq. (lA.ip . so we can neglect the time derivative in (1A.2|) . In 
addition, we assume that D, D c , a, and k are constants, and recast all variables in 
dimensionless form: t — > t t, r — > £q l D l l 2 r , p — > (D c /t ak)p, and c — > (D/k)c, where to 
is a typical timescale of the dynamics of p in Eq. flA.lj) . The resulting system is called 
the reduced Keller-Segel equations (00). 



Appendix B. Calculation of scalar products through Meijer G-function and 
T-function 

Calculation of scalar products in Section [5] requires to evaluate the integrals of the 
following type 



e-^y n+1 [ln(l + y 2 )]' 

(T+7? 



-dy, n, m, I e N, n > 0, 1 > 2, m > 0, (B.l) 



which by the change of variable x := 1 + y 2 and differentiation over the parameter a 
reduces to the following expression 
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And more generally, for n = and I >2: 
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A particular case m = is especially easy because G-function from fIB.ip reduces 

oo 

to the incomplete Gamma function T(s, z) = f t s ~ x e~ l dt as follows 
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A Taylor series expansion of ( IB. 51) for a — » gives for n = the following expressions 
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and the case n > is obtained by the differentiation of these expressions according to 
A Taylor series expansion of ( IB. 21) for a — > gives for n = the following expressions 
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The case n > is obtained by the differentiation of these expressions according to (IB.2j) . 
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